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A  model  is  presented  which  will  describe  the  time  evolution  of  oceanic  aerosol.  The  model  is 
designed  to  operate  on  a  desk  top  computer  (Tektronix  4052)  and  is  written  in  the  easy  to  under¬ 
stand  BASIC  language.  The  model  is  based  on  the  physical  processes  operating  on  a  column  of  air 
as  it  travels  with  the  wind  over  the  ocean.  Inputs  to  the  model  are  the  time  history  (or  forecasts) 
of  six  surface  observable  parameters:  air  temperature,  cloud  cover,  dewpoint,  windspeed,  sea  sur¬ 
face  temperature,  and  the  inversion  height.  The  output  of  the  model  is  a  series  of  aerosol 
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concentration  profiles  every  10  minutes  of  one  particular  dry  si/e  of  aerosol.  Provisions  arc  made 
for  tandem  runs  for  different  dry  sized  particles  so  that  size  distributions  of  aerosols  can  be  plotted 
as  they  evolve  over  time. 

The  model  takes  into  account  the  production  of  oceanic  aerosol  at  the  sea  surface  by  white 
caps  and  spray  and  the  changes  in  droplet  sizes  as  a  function  of  the  relative  humidity  profile  as 


well  as  the  mixing  processes  as  it  changes  with  both  time  and  altitude. 
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A  TIME-DEPENDENT  OCEANIC  AEROSOL  PROFILE  MODEL 


INTRODUCTION 

Marine  aerosols  have  always  been  of  concern  to  man.  In  the  early  days  of  seafaring  by  man.  opti¬ 
cal  visibility  was  of  major  concern  for  navigation  in  commerce  and  warfare.  Visual  observations  are 
important  to  make  navigational  fixes  on  shoreline  landmarks  as  well  as  for  sighting  the  stars.  The 
effects  of  visibility  during  wartime  played  an  important  role  in  the  outcome  of  many  naval  engage¬ 
ments.  The  advent  of  radar  in  World  War  II,  however,  did  away  with  this  strong  dependence  on  visual 
observations  because  radar  band  wavelengths  are  not  affected  by  the  existence  of  micron  sized  aerosol. 
With  the  use  of  radar  the  nav  igator  can  scan  the  distant  shoreline  for  known  land  marks.  Radio  naviga¬ 
tion  can  provide  fixes  to  unheard  of  precision  even  in  the  open  ocean.  Communication  with  other 
ships  and  with  land  bases  became  routine  with  radio  waves.  All  of  these  twentieth  century  advances 
using  radio  waves  are  essentially  independent  of  the  aerosol  content  of  the  marine  atmosphere. 

The  sophistication  of  the  modern  age,  however,  has  come  full  circle.  There  is  now  more  depen¬ 
dence  on  the  aerosol  content  of  the  atmosphere  because  different  wavelengths  of  electromagnetic  radia¬ 
tion  are  being  exploited  by  the  current  weapons  designers  due  in  part  to  the  technological  explosion  ol 
ideas  that  has  occurred  during  the  last  several  decades.  Thus  when  infrared  and  optical  radiation  are 
involved,  a  whole  area  of  its  interaction  with  atmospheric  aerosol  again  becomes  important.  This  is  par¬ 
ticularly  interesting  because  transmission  of  energy  through  an  aerosol  cloud  is  a  function  ol 
wavelength. 

The  Navy  has  a  need  to  predict  what  the  effect  of  meteorology  will  be  on  electro-optical  transmis¬ 
sion  characteristics  of  the  marine  atmosphere  and  shoreline  atmospheres  anywhere  throughout  the 
world  and  at  any  lime.  Obviously  this  problem  can  be  solved  only  if  a  suitable  maritime  aerosol  model 
is  developed  which  can  utilize  for  inputs  the  synoptic  scale  weather  information  that  is  gathered  rou¬ 
tinely  and  synthesized  by  various  large-scale  numerical  weather  predictive  systems  and  by  direct  satellite 
observations  of  the  globe  with  sensors  of  various  wavelength  emissions  and  the  appropriate  processing 
of  these  data. 

It  is  an  observation  of  seafarers  that  conditions  do  not  slay  constant  for  long  periods  of  lime. 
There  are  periods  of  calm,  and  then  come  the  storms  In  addition  the  environmental  history  of  an  air 
parcel  which  is  observed  on  a  ship  at  any  point  in  the  world's  ocean,  most  probably  is  vastly  different 
from  the  environmental  history  of  the  observing  platform  itself. 

Observational  data  of  aerosols  at  sea  [l!  show  that  at  intermediate  and  large  aerosol  sizes,  the 
aerosol  concentrations  increase  markedly  during  high  wind,  high  sea  slate  periods  and  lower  during 
calm  periods.  It  is  during  these  high  sea  state  periods  that  the  oceanic  component  of  the  aerosol  size 
spectrum  is  most  noticeable. 

The  purpose  of  this  report  is  to  describe  a  lime-dependent  aerosol  model  which  will  predict  the 
behavior  of  aerosol  throughout  the  planetary  boundary  layer  in  both  calm  and  white  water  conditions 
and  in  the  transitions  between  them.  This  model  will  respond  to  changes  in  the  marine  inversion 
height,  stability,  sea  state,  relative  humidity  profile,  and  wind  speed. 
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PHYSICAL  BASIS  OF  THE  MODEL 


These  processes  are  the  result  of  whiteeaps  and  spray  produced  by  the  action  of  wind  on  the  sea. 
Figure  1(a)  shows  a  cumulative  probability  curve  of  wind  speed  for  both  summer  and  winter  in  the 
North  Atlantic.  These  curves  indicate  that  90%  of  the  time  the  minimum  wind  speed  required  to  pro¬ 
duce  whiteeaps  is  exceeded  in  the  ocean.  Thus  droplets  from  white  water  could  be  a  very  important 
part  of  the  marine  aerosol  in  this  part  of  the  world. 


Central  to  the  ultimate  aims  of  this  effort  is  the  desired  ability  to  be  able  to  predict  the  future 
characteristics  of  the  aerosol  size  distribution  as  it  will  be  affected  by  generation  rates,  fallout,  and 
diffusion.  A  time-dependent  model  of  the  evolution  in  concentration  at  all  levels  of  the  boundary  layer 
of  a  single  dry  size  class  of  aerosol  is  the  means  by  which  this  need  will  be  met.  The  model  will 
describe  the  marine  aerosol  produced  by  processes  at  the  air-sea  interface.  The  droplets  are  composed 
of  water  and  sea  salt  and  are  therefore  very  hygroscopic  and  will  change  sizes  as  the  ambient  relative 
humidity  changes.  It  is  also  well  known  that  in  the  typical  marine  atmosphere  the  relative  humidity  is 
about  80%.  There  are  certainly  many  times  when  the  relative  humidity  is  about  80%  and  even  above 
90%.  Consider  for  example  the  plot  shown  in  Fig.  1(b).  Here  is  plotted  the  cumulative  probability 
curve  of  the  measured  relative  humidity  for  areas  such  as  the  North  Atlantic.  The  model  must  also 
take  these  changes  in  size  into  account. 


The  simples*  case  of  the  partial  differential  equations  which  will  adequately  describe  our  processes 
on  particles  of  dry  size  class  i  is  shown  in  Eq.  (1). 


dN, 

dt 


IkU) 


a/v, 

dz 


(1) 


This  equation  describes  the  processes  acting  on  a  column  of  air  moving  with  a  constant  velocity  in  the 
direction  of  the  wind.  This  equation  describes  the  time  history  of  the  concentration  of  particles  of  the 
size  class  /  at  an  altitude  c  above  the  sea  surface.  There  are  three  physical  processes  assumed  acting  on 
the  particles  and  they  are: 


(1)  generation  rate  at  the  sea  surface 

(2)  eddy  diffusion 

(3)  gravitational  fall. 


The  generation  rate  will  be  reflected  by  the  lower  boundary  condition.  The  mean  micrometeorological 
state  of  the  atmosphere  is  described  by  an  altitude-dependent  eddy  diffusion  coefficient  represented  by 
K  (z).  The  fall  rate  of  the  particles  in  class  /  is  represented  by  V,  and  will  be  a  function  of  the  ambient 
relative  humidity  which  will  in  general  be  different  for  each  altitude. 

To  adequately  represent  the  time-dependent  solution  of  this  differential  equation,  certain  boun¬ 
dary  conditions  must  be  defined.  At  time  equals  zero,  the  initial  concentrations  of  particles  throughout 
the  column  must  be  specified.  If  we  are  interested  in  only  observing  that  portion  of  the  total  aerosol 
spectrum  which  is  being  locally  generated  by  wave  action,  then  we  may  specify  t Ve  initial  conditions 
such  that  there  is  a  negligible  amount  of  aerosol  throughout  the  column  at  the  slat*  of  the  calculation 
time.  This  is  equivalent  to  looking  at  the  cases  where  the  continental  component  of  the  atmospheric 
aerosol  is  always  equal  to  zero.  In  addition,  surface  boundary  conditions  must  be  specified.  In  this 
simple  one-dimensional  case,  only  the  conditions  at  the  top  and  the  bonom  of  our  column  of  air  must 
be  specified.  The  upper  boundary  condition  is  easi’y  taken  care  of  if  the  height  of  our  air  column  is 
substantially  greater  than  the  height  of  the  inversion  layer  specified  in  the  K(z)  function.  Since  we  are 
only  looking  at  sources  of  aerosol  at  the  sea  surface  and  the  capping  inversion  specifies  no  diffusion  of 
panicles  through  the  inversion,  then  there  is  no  problem  in  arbitrarily  keeping  the  upper  boundary  con¬ 
stant  in  value  as  time  progresses.  On  the  other  hand,  the  white  wat.er  processes  at  (he  sea's  surface  will 
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greatly  affect  the  concentrations  at  the  lowest  levels  of  the  column  and  in  addition  may  even  vary  as  a 
function  of  time.  To  simulate  this  process,  a  flux  boundary  condition  is  maintained  at  the  sea  surface 
which  will  be  varied  as  a  function  of  the  meteorological  input  parameters. 


EDDY  DIFFUSION 

The  problem  of  defining  the  eddy  diffusion  coefficient  as  a  function  of  altitude  is  taken  care  of  in 
this  model  by  using  the  work  of  O'Brien  [21,  where  he  has  constructed  a  functional  form  for  K(z) 
which  is  continuous  in  both  its  value  and  its  derivative.  The  general  characteristics  of  this  function  can 
be  matched  to  the  atmospheric  conditions  by  parameter  setting.  Of  particular  importance  in  this  regards 
is  the  parameter  describing  the  height  of  the  inversion  over  a  typical  open  oceanic  area.  This  altitude 
will  correspond  to  the  top  of  the  marine  layer.  This  is  the  altitude  at  which  we  wish  to  make  the  eddy 
diffusion  function  to  become  small  or  zero  under  neutral  and  unstable  conditions. 

The  process  this  model  goes  through  to  determine  the  profile  parameters  necessary  to  define  the 
functional  form  of  the  eddy  diffusion  coefficient  starts  with  the  calculation  of  various  micrometeorologi- 
cal  scalar  parameters  based  on  the  six  bulk  meteorological  input  parameters.  These  parameters  are 
stored  in  a  storage  array  for  future  reference.  This  process  commences  in  task  2  of  the  program  with 
the  calculation  of  the  drag  coefficient,  CIO.  This  parameter  (stored  in  P(6)  of  the  scalar  storage  array) 
is  linearly  related  to  the  measured  wind  velocity  in  the  functional  form  suggested  by  Deacon  [3,4]. 

C, o-  0.0011  +  4x  10” 5  x  Ux o,  (2) 

where  U j0  is  the  wind  speed  at  10  m  measured  in  meters  per  second. 

The  friction  velocity,  U *  is  next  calculated  from  the  input  wind  speed  at  10  m  and  the  previously 
calculated  drag  coefficient.  This  parameter  is  stored  in  P(8)  with  the  units  of  meters  per  second. 

u.  =  Vc7o  x  £/|0-  (3) 

The  dynamic  roughness  of  the  sea  surface,  Z0,  is  next  calculated  from  the  friction  velocity  in  the 
manner  of  Charnock  [5]. 

Z0  =  C, o  x  Ufa  (4) 

This  parameter  is  stored  in  P(9)  with  the  units  of  meters. 

Next  the  bulk  derived  Richardson's  number  is  calculated  from  the  gradient  in  virtual  potential 
temperature  and  wind  speed  at  the  associated  10-m  level  [6]  by: 

Rivm  =  3.55(0v,o-  0vo)/t/,2o  (5) 

where  the  influence  of  the  water  vapor  gradient  on  stability  is  accounted  for  by  the  gradients  in  the  vir¬ 
tual  potential  temperature  calculated  from  the  sea  surface  temperature,  the  air  temperature,  and  the 
dewpoint  in  the  standard  manner.  Once  the  bulk  Richardson’s  number  is  calculated,  it  is  stored  in  the 
location  P(10).  This  number  will  be  used  in  later  determinations  to  describe  the  eddy  diffusion 
coefficient  profile. 

The  actual  stability  of  the  atmospheric  boundary  layer  will  be  determined  by  the  ratio  of  the  refer¬ 
ence  altitude  Z|0  to  the  Monin-Obukhov  length,  L.  Negative  values  of  this  ratio  rl0/L  indicate  instabil¬ 
ity  whereas  positive  values  indicate  stability.  The  dimensionless  stability  parameter  z{0/L  is  calculated 
by  the  method  of  Barker  and  Baxter  [7]  from  the  bulk  Richardson's  number  and  is  stored  for  later 
reference  in  P(7). 


\  K I  RU'UK  I  #\*h 

The  calculation  of  ihe  O'Brien  1 2 1  polynomial  form  of  ihc  eddy  diffusion  coefficient.  Ate), 
requires  values  ol  the  eddv  difiusion  at  1  m.  A (1  I  and  its  derivative.  A(l)  at  the  l-m  level  a  t  II  and 
AO)  are  computed  from  the  formula: 

A  iz  )  -  0.35  U.  z/4>iz/L  ).  (6) 

were  <Mr/A)  is  the  semiempirical  stability  function  for  heat  transport  |8l 

For  unstable  conditions  the  depth  of  the  planetary  boundary  layer  is  taken  as  the  height  of  the 
capping  inversion  stored  as  one  of  the  input  time  parameters.  During  stable  conditions,  Ciarke  (9] 
empirically  determined  that  turbulent  mixing  approaches  zero  at  a  level: 

//  »  0.23  x  UJ  {Corioles  parameter].  (7) 

For  very  stable  conditions  (Rib  >  1/4.7),  no  turbulent  mixing  is  predicted  by  the  model  and  Kiz)  is 
everywhere  negligible. 

RELATIVE  HUMIDITY  EFFECTS 

The  dynamic  aerosol  model  is  able  to  determine  the  effect  of  droplet  growth  in  near  saturation 
conditions.  Information  on  the  relative  humidity  profile  of  the  atmosphere  is  obtained  from  model 
predictions  based  on  surface  shipboard  observations  [10].  The  dynamic  aerosol  model  assumes  a  profile 
based  on  the  previously  declared  inversion  height  and  surface  meteorological  measurements. 

In  the  model  the  hygroscopic  particles  produced  at  the  sea  surface  are  always  in  size  equilibrium 
with  the  relative  humidity  of  their  surroundings.  When  the  computer  first  asks  what  size  particle  is 
desired  to  follow,  it  also  asks  at  what  "reference"  relative  humidity  is  to  be  used  Thus  with  the  formu* 
lations  of  Fitzgerald  [11]  it  determines  the  dry  size  of  the  particle  and  thus  can  predict  to  what  size  it 
would  grow  in  any  other  subsaturated  relative  humidity  environment  in  which  they  would  be  placed. 
Since  these  particles  are  all  of  the  same  dry  size  no  matter  how  large  they  grow,  we  can  at  any  altitude 
talk  about  the  number  of  particles  which  have  the  same  dry  size.  Thai  is  the  particles  are  always 
identifiable  by  their  dry  size.  Therefore  even  though  the  relative  humidity  changes  we  can  always  keep 
track  of  the  particular  dry  sized  particles.  Thus  given  any  size  distribution  of  dry  sized  particles  we  can 
easily  determine  the  size  distribution  which  this  population  would  have  if  it  were  immersed  in  other 
relative  humidity  environments. 

The  fall  velocity  of  the  droplets  is  the  only  quantity  in  which  the  defining  partial  differential  equa- 
lion  is  affected  by  the  droplet  size  and  therefore  most  strongly  by  the  relative  humidity  profile.  Thus 
this  is  accounted  for  in  the  model  by  adjusting  the  fall  velocity  for  each  different  level  to  that  which  a 
particle  of  the  initial  dry  size  would  have  if  it  were  grown  to  the  size  indicated  by  formula  for  the  rela¬ 
tive  humidity  at  that  level. 

To  implement  this  feature,  the  profile  of  relative  humidity  must  be  known.  Such  information  is 
not  ordinarily  known  unless  a  sounding  has  been  done.  Gathman  [10],  however,  has  devised  an  empiri¬ 
cal  model  for  predicting  the  relative  humidity  profile  given  only  the  sea  surface  meteorological  mea¬ 
surements.  The  main  features  of  this  model  are  incorporated  into  this  time-dependent  oceanic  aerosol 
model,  and  its  structure  is  based  on  the  hourly  surface  meteorological  inputs. 

The  sedimentation  rate  for  water  droplets  in  air  is  determined  by  the  Stokes-Cunningham  relation¬ 
ship  where  the  settling  velocity  expressed  in  meters  per  second  is  given  by: 

V  =  3004  x  D2  x  (I  +  1.6  xIO  5//>).  (8) 

where  D  is  the  diameter  of  the  particle  in  centimeters  at  the  ambient  relative  humidity  (12) 
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LOWER  BOUNDARY  CONDITION 

The  model  uses  as  ihe  lower  boundary  condition,  the  upward  flux  of  particles  which  occur  over 
salt  water  for  that  particular  sea  slate.  This  approach  essentially  uses  an  empirical  wind-dependent  rela¬ 
tionship  to  define  the  aerosol  flux  which  is  necessary  for  the  profile  calculation  of  the  oceanic  com¬ 
ponent  of  the  atmospheric  aerosol. 

The  works  of  Woodcock  [131  who  obtained  salt  particle  distributions  as  a  function  of  wind  speed 
and  of  Eriksson  [141,  who  computed  the  fall  velocities  for  salt  particles  in  equilibrium  with  a  91.4% 
relative  humidity  are  the  basis  by  which  Blanchard  [15]  obtained  a  set  of  curves  that  give  (under  the 
assumption  of  steady  state)  the  spectra  of  the  number  of  ocean  produced  drops  leaving  the  sea  surface 
per  unit  area  per  second  per  0.2xLog(D)  interval  (where  D  is  the  droplet  diameter).  Blanchard  pro¬ 
duced  different  spectra  for  different  wind  speeds.  This  family  of  curves  can  be  fitted  by  a  set  of  log¬ 
normal  distributions,  the  three  coefficients  of  which  are  wind  dependent  and  can  be  themselves 
represented  by  simple  functions  of  the  wind  speed.  The  fitted  spectra,  shown  in  Fig.  2,  closely  resem¬ 
ble  Blanchard's  figure  and  will  be  the  basis  of  this  model. 


Fig  2  —  The  aerosol  flux  size  distribution  plotted  as  a  log 
normal,  wind-dependent  function  fitted  to  Blanchard's  data 


Thus  by  knowing  the  wind  speed,  we  can  compute  the  three  coefficients  which  describe  the  paitic- 
ular  log-normal  size  distribution  of  the  droplet  flux.  By  knowing  the  size  our  particular  dry  sized  parti¬ 
cles  have  at  the  relative  humidity  of  91.4%  we  can  then  easily  obtain  the  surface  flux  boundary  condi¬ 
tion. 


The  extrapolation  of  these  curves  to  low  wind  speeds  shows  that  at  wind  speeds  below  4.5  knots, 
the  amplitude  parameter  of  the  log-normal  distribution  drops  to  zero.  This  result  is  to  be  expected 
because  of  the  lack  of  whitecaps  at  these  speeds.  The  curves  also  show  a  general  broadening  of  the  dis¬ 
tribution  with  increasing  wind  speed  as  well  as  the  shifting  of  the  maximum  concentration  size  to  larger 
values  at  the  higher  wind  speeds. 
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The  droplets  in  this  distribution  in  the  diameter  range  of  from  6  to  60  microns  are  presumably  the 
results  of  the  bursting  of  bubbles  unci  the  spray  from  waves  and  do  not  include  the  very  small  droplets 
which  result  from  the  fragmentation  of  the  bubble  film  in  the  bubble  bursting  process. 

MODEL  IMPLEMENTATION 

This  model  deals  with  the  expected  concentrations  of  oceanic  aerosol  of  a  particular  dry  si/e  of 
particle.  As  the  relative  humidity  cannot  generally  be  expected  to  be  invariant  with  altitude  we  are  then 
talking  about  the  concentrations  of  different  real  size  droplets  at  each  of  the  levels  of  our  model  grid 
because  the  relative  humidity  can  be  expected  to  be  different  at  each  of  these  levels.  The  fact  remains, 
however,  that  anv  of  our  droplets  if  dried  out  would  approach  a  common  size. 

To  tag  each  of  the  droplets  so  that  the  computer  can  always  keep  track  of  the  particular  droplet  we 
are  studying,  both  a  diameter  and  the  relative  humidity  must  be  known.  Thus  as  a  computer  input,  a 
diameter  and  a  relative  humidity  at  which  this  diameter  is  referenced  must  be  given.  Once  this  is 
known,  a  dry  size  is  calculated  and  all  processes  affecting  this  set  of  constant  dry  sized  particles  are 
computed. 

The  model  is  designed  to  predict  the  concentration  of  only  one  dry  size  of  particle  at  anv  particu¬ 
lar  run.  This  is  a  limitation  on  the  size  and  speed  of  the  computer  being  used.  Size  distributions  can 
be  obtained  by  running  the  model  several  times  with  the  same  set  of  input  parameters  except  that 
different  dry  sizes  of  particles  will  be  given  to  the  computer  to  be  considered.  As  the  concentrations  at 
all  levels  of  these  various  dry  sized  particles  are  predicted  by  the  model  as  well  as  their  ambient  relative 
humidities,  the  size  distribution  at  any  time  and  at  any  altitude  can  be  calculated  by  tandem  computer 
runs. 

If  a  log  normal  distribution  is  adequate  to  represent  the  size  distribution,  then  only  three  runs 
need  be  done  as  these  type  of  distributions  can  be  completely  defined  by  the  concentrations  of  three 
sizes  of  particles  (any  three  sizes  of  particles). 

A  rerun  provision  is  included  in  this  program  which  will  remember  all  of  the  meteorological  his¬ 
tory  vectors  in  order  to  provide  identical  runs  for  different  sized  droplets.  A  range  of  droplet  sizes 
which  this  models'  input  data  is  matched  is  from  a  2  to  90  microns  diameter,  at  relative  humidities  of 
approximately  90%.  Sources  of  droplets  from  the  ocean  surface  such  as  film  droplets,  jet  droplets,  and 
spray  fall  with  in  this  range.  Thus  asking  the  computer  to  look  at  droplets  outside  of  this  range  is 
prohibited  by  the  software. 

The  model  results  depend  on  certain  meteorological  measurements  as  inputs  and  boundary  condi¬ 
tions  which  are  required  information  in  order  for  the  model  to  operate.  Hourly  meteorological  measure¬ 
ments  are  assumed.  When  used  to  look  backwards  at  the  atmospheric  aerosol  history  actual  hourly 
meteorological  records  can  be  given  to  the  computer.  The  computer  assumes  the  "current” 
meteorological  environment  is  a  good  approximation  until  told  of  future  changes  during  the  input  pro¬ 
cess.  Standard  weather  forecasts  or  the  use  of  more  sophisticated  model  forecasts  may  also  be  used  to 
provide  the  meteorological  history  parameters  for  forward-looking  predictive  runs. 

Various  calculations  based  on  the  measured  meteorological  parameters  are  used  to  prov  ide  those 
specifically  needed  by  the  model  and  are  included  as  subroutines  which  can  be  changed  as  necessary. 
For  instance  the  mierometeorological  estimates  used  in  the  calculation  K  iz)  based  on  O'Brien  [2]  can 
be  easily  updated  as  research  in  that  field  progresses. 

The  inversion  height  is  considered  a  meteorological  parameter  and  entered  by  hand.  These  may 
come  from  several  sources  such  as  an  acoustic  sounder,  a  lidar,  a  radiosonde,  a  model,  or  from  a  clima¬ 
tology 
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The  mode  I  in  implemented  in  the  BASK'  language  on  the  lektromx  4052  desk  top  graphic  com¬ 
puter  equipped  with  disc  storage.  1  he  model  lias  certain  characteristics  ot  i tic  machine,  but  these 
features  are  not  absolutely  necessary  to  the  basic  structure  and  How  of  the  model  I  he  program, 
without  the  added  special  features,  should  be  runable  with  few  modifications  on  am  BASK  language 
computing  machine  with  adequate  memory 

The  speed  of  computations  will  of  course  \ary  from  machine  to  machine  and  could  make  main 
small-scale  computer  systems  require  an  excessive  amount  of  time  for  computation  On  the  other  hand 
if  speed  of  computation  is  of  utmost  importance,  the  program  could  be  converted  to  a  compiler-like 
language  and  run  with  small  modifications  on  a  last  large-scale  computer. 

The  method  used  in  the  solution  of  this  initial  value  problem  with  one  spatial  dimension  is 
described  in  Young  and  Gregory  1 16]  as  the  forward  difference  method.  To  increase  the  grid  spacing 
close  to  the  sea's  surface  where  the  largest  spatial  variations  in  the  various  quantities  takes  place,  a  vari¬ 
able  grid  is  introduced  for  24  levels  used  to  describe  this  model. 

The  computational  stability  of  this  method  of  solution  depends  on  the  mesh  ratio  r  which  K 
defined  as  the  ratio  of  the  time  step  to  the  square  of  spatial  grid  spacing.  It  can  be  shown  that 
A  (r)  x  r  must  be  less  than  or  equal  to  1/2  in  order  to  assume  satisfactory  results  in  the  computation. 
In  the  variable  grid  scheme  used  here,  the  lowest  grid  spacing  is  usually  the  critical  one  which  limits  the 
speed  of  computation.  The  program  is  set  up  so  that  the  Tektronies  4052  will  normally  process  24 
hours  of  computational  time  in  1  hour  of  wall  clock  time.  Given  a  particular  machine  the  processing 
speed  can  be  improved  only  by  degrading  the  grid  spacing  such  that  the  stability  conditions  in  the  mesh 
ratio  remains  valid.  However  to  maximize  the  speed  of  computation,  a  new  maximum  time  step  inter¬ 
val  is  calculated  after  each  new  K  (:)  is  determined. 

MODEL  INPUTS 

The  interacting  capabilities  of  the  4052  computer  are  utilized  to  allow  the  computer  to  ask  ques¬ 
tions  of  the  operator  so  that  all  of  the  required  input  parameters  are  satisfied.  The  program's  meteoro¬ 
logical  input  data  are  stored  in  the  form  of  time  lists,  with  estimated  or  measured  values  for  each  of 
these  meteorological  quantities  required  for  each  "model  time"  hour  that  the  calculations  are  made.  In 
inputting  this  data,  the  concept  of  persistence  is  used  in  that  an  initial  value  is  required  for  the  first 
model  time  hour  of  each  parameter  and  these  initial  values  are  assumed  for  the  remaining  23  hours 
unless  they  are  specifically  changed  for  some  later  times.  If  they  are  changed,  the  new  values  are  used 
until  the  twenty-fourth  hour  unless  it  is  again  changed.  There  are  six  meteorological  time  lists  used  in 
the  model,  and  they  are  listed  below  together  with  their  required  units  for  the  input  process. 

(1)  ,4  1  =  Air  temperature  history  (°C) 

(2)  Cl  =  Cloud  cover  history  (tenths) 

(3)  D  \  =  Dew  point  history  (°C) 

(4)  H 9  =  Inversion  base  height  history  (m) 

(5)  V 9=  Wind  speed  history  (m/s) 

(6)  .5*0=  Sea  surface  temperature  history  (°C) 

In  addition  to  the  time  lists,  certain  parameters  concerning  the  calculation  are  required  before  the 
program  can  begin.  The  general  geographical  latitude  is  required  (in  degrees)  as  well  as  the  droplet  size 
of  interest  (expressed  as  a  diameter  in  microns)  together  with  a  reference  relative  humidity  at  which 
the  desired  droplet  size  is  associated.  This  latter  information  is  used  to  determine  the  actual  dry  size  of 
the  particle  so  that  its  changes  in  size  with  respect  to  the  various  relative  humidities  encountered  in  the 
mode!  can  be  accurately  determined. 
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I  he  mode!  nukes  a  numhu  ol  t alculaiions  based  mi  the  mcleMit.louie.il  input  cl.i  1.1  Pioliles  ol 
the  cdd>  ilitUiMon  coefficient  and  the  relative  humulitv  ,i>  well  ,i'  the  value  «•!  the  surlace  Hus  .me  based 
emuelv  mi  the  lime  list  data  and  arc  calculated  imtiailv  and  are  updated  each  "model  time'  hour  there¬ 
after  in  provide  the  eh  immic  framework  mi  which  the  phvMcai  processes  of  faf/ou;  and  edd>  dittusum 
o per  ate 

l  he  determination  of  these  three  ciuaniities  contains  much  ot  the  ph>siws  mi  which  the  model 
depends  and  in  turn  alfeets  the  si/e  of  the  droplets  and  thus  the  falling  velocities,  the  abiJih  of  the 
houndarv  laver  to  redistribute  the  droplets  which  are  in  it.  and  the  oceanic  source  o|  these  droplets 
from  white  water  phenomenon  The  methods  for  determining  these  quantities  tmm  the  surl.ue 
meteorological  inputs  are  based  on  other  works  l>resumabl>  as  research  m  these  areas  continues,  an 
improvement  in  this  model  can  be  obtained  bv  simply  upgrading  the  various  subroutines  used  in  the 
model  These  data  are  stored  as  elements  of  parameter  vector  l\  and  the  definitions  ot  these  elements 
are  shown  in  Table  1 

Table  I  —  Definition  of  the  Storage  Array  P 


P(D  .  DROPLET  DIAMETER  (MICRONS) 

PI  2)  -----  REFERENCE  RELATIVE  HUMIDITY  Of  SI/E  ME  AS  \  ) 

Pi 3)  .  DROPLET  DRY  DIAMETER  (MICRONS) 

PC  5 )  .  LATITUDE 

P(6)  .  DRAG  COEFFICIENT,  CIO 

PI  7 )  .  Z/L 

P(8 >  .  FRICTION  VELOCITY,  U* 

Pi 9)  .  SURFACE  ROUGHNESS,  70 

PI  10)  .  BULK  RICHARDSON'S  NUMBER 

PHD  -----  FRICTION  POTENTIAL  TEMPERATURE  Till  1  A“ 

Pi  12)  -----  FRICTION  MIXING  RATIO.  R*  OR  Q' 

Pi  13)  .  HEIGHT  OF  SUPER  ADIABATIC  LAYER.  ZSA 

PI  14)  .  HEIGHT  OF  CONDENSATION  LAYER.  /CON 

P(  15)  .  F  LUX  OF  DROPLETS  F  ROM  THE  SURFACE 

PI  16)  .  HOUR  COUNTER 

P(  1 7 )  .  PARAMETER  UPDATE  TRIGGER  (NORMALLY  -0) 

PI  18)  .  DATA  STORAGE  BLOCK  LOCATOR 

PI  19)  .  DELTA  TIME  STEP 

P( 20)  TRIGGER  VALUE  FOR  DISPLAY  AND  OPERATOR  DIRECTION 

P ( 2 1 )  TRIGGER  FOR  AUTOMATIC  DATA  STORAGE 

P(22)  .  TIME  HOURS 

PI 23)  TIME  #  OF  10  MINUTE  INTERVALS 

PI 24)  .  TIME  SEC 


FLOW  CHART 

The  flow  of  the  program  is  illustrated  in  Fig.  3  where  the  main  features  of  the  model  steps  are 
shown  as  blocks.  Because  certain  of  the  steps  arc  lengthy  and  require  much  memory  space,  the  pro¬ 
gram  is  broken  down  into  tasks  for  convenience.  The  programs  for  these  tasks  and  their  required  sub¬ 
routines  are  stored  on  disc  and  are  copied  into  an  overlay  section  of  the  main  program  when  they  are 
needed  The  portions  of  the  program  flow  chart  which  are  identified  with  the  various  tasks  are  outlined 
in  the  figure.  Most  of  the  program  calculation  lime  is  spent  in  task  5  so  that  the  overlaying  does  not 
take  up  too  much  time. 
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The  main  product  of  the  program  is  a  disk  file  written  as  a  random  access  file  which  stores  the 
complete  aerosol  profile  every  10  minutes  of  "model  time."  The  information  in  this  file  can  he  eusilv 
accessed  by  other  analysis  programs  to  describe  the  evolution  of  the  aerosol  concentration  of  the  partic¬ 
ular  dry  si/es  asked  for  at  the  beginning  of  the  calculation. 

The  basic  model  describing  mixing  by  eddy  diffusion  and  the  gravitational  fall  of  droplets  is 
affected  by  the  meteorological  measurements  used  in  the  functional  form  of  the  profile  of  the  eddy 
diffusion  coefficient  A  (r)  and  the  fall  velocity  V iz)  of  the  droplets  themselves,  which  depends  on  the 
size  of  the  droplets,  which  in  turn  depends  on  the  ambient  relative  humidity  of  the  atmosphere  The 
lower  boundary  condition  in  the  model  presented  here  V  1(1)  is  also  related  empirically  to  wind  speed 

Because  the  predictions  of  the  model  do  not  in  any  way  affect  the  meteorological  values  of  the  air 
column  being  modeled,  there  is  no  feedback  process  occurring.  Hourly  observations  are  usually  ade¬ 
quate  to  describe  the  meteorological  features  of  mesoscale  processes,  and  this  frequency  of  updating  is 
assumed  adequate  for  the  time-dependent  oceanic  aerosol  profile  model. 

Ever>  model  hour  the  meteorologically  sensitive  profiles  A'(r),  V ( r ) .  and  A'l(r)  are  updated 
from  the  best  estimates  of  the  basic  meteorological  input  parameters  while  the  calculation  continues.  In 
this  way  the  lengthy  calculations  of  these  values  are  not  repeated  more  often  than  necessary. 

MODEL  PERFORMANCE 

The  complete  testing  of  a  comprehensive  numerical  geophysical  model  is  an  extremly  difficult  pro¬ 
cess  because  of  the  large  number  of  combinations  of  parameters  which  are  allowed  as  inputs.  \  given 
field  test  may  be  a  basis  for  testing  the  model  for  a  finite  set  of  input  conditions,  but  there  is  alwa> s  a 
chance  that  one  certain  set  of  parameters  will  cause  completely  erroneous  results.  Also,  the  numerical 
solution  of  a  problem  causes  further  uncertainties  when  finite  grid  spacings  and  lime  steps  are 
employed. 

One  approach  of  checking  for  the  errors  introduced  into  the  numerical  solution  is  to  find  a  case 
that  can  be  solved  by  analytical  means  and  then  to  check  the  numerical  solution  with  the  analvtical 
solution  to  see  how  valid  the  numerical  approximation  is. 

Other  tests  that  can  be  done  to  strengthen  the  plausibility  of  the  model's  results  are  to  run  very 
specific  cases  where  the  outcome  is  generally  known.  All  inputs  in  these  cases  can  be  simplified  and 
made  to  remain  constant  except  the  particular  parameter  for  which  the  model's  sensitivity  is  desired  to 
be  known. 

The  defining  partial  differential  equation  (Eq.  (1))  may  be  solved  analytically  for  steady  state 
cases.  The  numerical  solution  of  the  time-dependent  oceanic  aerosol  model  approaches  the  steady  state 
solution  (if  the  boundary  conditions  remain  constant)  as  time  approaches  very  large  numbers.  Thus  a 
check  of  the  digital  technique  can  be  obtained  by  comparisons  with  a  steady  state  analytic  solution. 

For  this  purpose  the  subroutines  defining  the  profiles  of  eddie  diffusion  coefficient  and  the  relative 
humidity  were  changed  in  the  numerical  model  to  correspond  with  integratable  functions  for  the  ana¬ 
lytic  solution.  To  this  end  the  relative  humidity  profile  is  also  held  constant  at  the  reference  value  so 
that  the  fall  velocity  would  not  be  a  function  of  altitude  but  only  of  the  droplet  dry  size. 

The  altitude  dependence  of  the  proposed  eddie  diffusi  n  coefficient  for  the  comparison  between 
the  digital  computer  solution'and  the  analytical  model  was  chosen  to  be: 

K(z)  =  K  x  :  +  g.  (9) 


where  K  and  x  are  arbitrary  constants. 
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In  the  steady  state  case  dX/dt= 0  and  therefore  Eq  (!)  becomes 

K  (z)  x  dS/dz  +  v  x  \  =  constant.  HO) 

If  AH)  is  defined  as  in  Eq.  (9)  and  v  is  held  constant,  the  equation  integrates  to: 

\  (z )  -  CM  x  (A  x  :  +  K)"  (2  HI) 

where  Cl  and  C 2  are  constants  of  the  integration. 

In  the  comparison  between  the  analytical  and  the  numerical  solution,  we  defined  A  =  2  and 
g  =  0.01,  and  solved  the  problem  for  a  40  micron  droplet  in  an  ambient  80%  relative  humidity  This 
droplet  has  a  fall  velocity  of  0.0482  m/s. 

When  a  computer  run  of  five  and  a  half  model  hours  was  done,  the  digital  solution  had  reached  a 
steady  state  condition  and  the  integration  constants  Cl  and  C2  were  evaluated  from  the  numerical 
model  levels  of  1 1  and  594  m.  Figure  4  is  a  plot  of  the  error  in  percent  between  the  analytical  solution 

and  the  numerical  s  * tion  plotted  as  a  function  of  altitude.  This  curve  shows  the  effects  of  large  grid 

spacing^  at  the  higher  altitudes  on  the  accuracy  of  the  solutions.  Likewise  the  larger  errors  at  the 
lowest  levels  probably  reflect  the  matching  of  boundary  conditions  at  the  boundary  itself. 


Fig  4  —  Error  between  an  analytical  solution  and  a  steady 
state  digital  solution  for  a  simplified  profile  of  relative 
humidity  and  eddy  diffusion 


A  second  example  in  which  the  performance  and  importance  of  the  time  dependent  oceanic  aero¬ 
sol  model  can  be  seen  is  the  response  of  the  atmosphere  to  a  sudden  increase  in  wind  speed  for  a  short 
duration.  In  the  carrying  out  of  this  test  the  subroutines  of  the  model  are  all  restored  to  their  normal 
form.  The  meteorological  parameters  are  set  to  be  thermally  stable,  and  the  relative  humidity  was 
forced  to  be  very  high  so  as  to  simulate  an  encounter  in  a  slightly  subsaturated  condition.  The  wind 
speed  input  was  set  to  zero  for  the  first  hour,  15  m/s  during  the  second  hour,  and  again  zero  for  the 
remainder  of  the  calculation.  In  this  case  of  course,  there  is  no  input  flux  unless  the  wind  speed  is 
high  Likewise  the  eddie  diffusion  parameter  is  very  low  except  for  the  mechanical  mixing  produced 
during  the  period  of  high  wind  speed.  Four  runs  were  made  of  this  case  for  the  sizes  of  20,  10,  5,  and 
2.5  microns  diameter  with  a  relative  humidity  of  90%. 
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The  results  of  the  first  500  min  of  this  run  are  plotted  in  Fig. 5.  These  three  dimensional  plots 
show  the  surface  defined  by  the  altitude,  concentration,  and  time  of  a  certain  dry  si/ed  aerosol  com¬ 
puted  by  these  runs.  The  vertical  axis  is  altitude  plotted  linearly.  The  horizontal  axis  is  the  linear 
measure  of  aerosol  concentration.  The  time  axis  goes  back  into  the  page  at  an  angle.  The  altitude-time 
plane  is  a  background  concentration  profile  of  \(:)  =  10  >  per  cubic  meter  which  is  represented  by  the 
vertical  lines  describing  constant  V(r)  profiles  every  10  min  and  the  slanted  lines  depicting  levels  of 
constant  altitude  which  show  the  variable  grid  levels  used  in  the  model.  For  comparison  purposes,  the 
scale  of  the  concentration  is  not  constant  but  is  adjusted  so  that  the  maximum  concentration  distance  of 
the  whole  surface  is  the  same  for  each  of  the  aerosol  size  classes.  The  absolute  values  of  these  concen¬ 
trations  vary  by  a  considerable  amount. 


Fig  5  —  Three-dimensional  surfaces  describing  ihe  time  and  space  evolution  of  various  size  droplets 


The  interesting  feature  of  these  surfaces  is  the  time  responses  of  the  atmosphere  to  each  of  the 
aerosol  sizes.  The  largest  aerosol  size  class  mirrors  the  wind  speed  function  and  does  not  affect  the 
atmosphere  above  50  m  while  the  smallest  size  class  appears  to  just  slay  in  the  atmosphere  once  intro¬ 
duced  into  it.  It  of  course  eventually  falls  out  but  at  a  much  slower  rate  than  the  largest  aerosol.  This 
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of  course  will  be  exhibited  in  a  changing  size  distribution  which  will  narrow  and  shift  to  smaller  sizes  as 
time  goes  on  in  the  event  there  are  no  generation  processes. 

This  shifting  of  the  size  distribution  can  be  seen  in  Fig.  6  where  the  atmosphere  was  filled  initially 
everywhere  with  the  same  log-normal  distribution  of  aerosol  shown  by  the  solid  line  in  the  figure.  The 
model  was  allowed  to  operate  for  1 1  hours  of  model  time  under  idealized  conditions.  These  conditions 
included  a  constant  relative  humidity  profile,  and  a  wind  speed  such  that  no  flux  was  introduced  from 
the  surface  but  yet  there  was  enough  eddy  diffusion  to  keep  the  boundary  layer  stirred.  The  resulting 
size  distribution  is  shown  by  the  crosses  in  the  figure.  This  calculation  shows  that  indeed  the  size  dis 
tribution  does  change  with  time  as  expected. 


DIAMETER  (microns) 

Fig  6  —  The  modification  of  an  aerosol  size  distribution  in 
a  nonequilibrium  slate.  The  curve  is  the  initial  size  distri¬ 
bution  throughout  an  atmospheric  column  The  crosses  are 
the  model  calculations  at  an  altitude  of  H  m  of  the  modified 
size  distribution  after  a  period  of  II  h  of  mixing  produced 
at  a  level  just  below  the  white  cap  threshold 


CONCLUSIONS 

The  performance  of  the  model  gives  plausible  results  when  certain  "setup"  cases  are  introduced  to 
it.  There  are  several  sets  of  field  data  available  which  provide  enough  of  both  the  aerosol  size  distribu¬ 
tions  and  the  necessary  meteorology  input  data  so  that  it  can  be  tested  on  a  case  study  basis.  These 
lengthy  case  studies  are  beyond  the  scope  of  this  report  but  are  presently  being  processed  and  prepared 
in  the  form  of  a  companion  report. 

The  convenience  of  operating  the  model  on  a  small  desk  top  computer  is  very  great.  Numerical 
experiments  can  be  carried  out  easily  and  inexpensively  and  the  results  plotted  with  ease.  As  the  origi¬ 
nal  desire  was  expressed,  the  most  productive  use  of  the  model  is  expected  to  be  the  ability  to  use  fore¬ 
casted  meteorological  data  as  well  as  historical  meteorological  data  to  determine  the  future  of  an  evolv  ¬ 
ing  vertical  profile  of  the  oceanically  produced  aerosol  size  distribution. 
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1  GO  TO  lOOOO 
B  REM  PRINT  MODEL  TIME 

9  P ( 20 ) -2 

10  RETURN 

12  REM  PLOT  N 1 

13  P ( 20 ) =3 

14  RETURN 

16  REM  PRINT  N1 

17  P<20>=4 

18  RETURN 

20  REM  SAVE  WORLD  ROUTINE 

21  P ( 20 ) =5 

22  RETURN 

10000  GO  TO  30090 


20000  REM - OVL  —  AREA 

20010  REM 
30000  REM 

30010  REM  ! - 

30020  REM  i  MAIN  PROGRAM 

30030  REM  I  MAIN 

30040  REM  I - 

30050  REM 

30055  REM - LOAD  OVL - AREA 

30060  GOSUB  30220 
30070  GO  TO  20060 

30080  REM - MAIN  PROGRAM 


30090  INIT 
30100  Z$="" 

30110  U=1 

30120  REM  - 

30130  REM  ; 

30140  REM  ! 

30150  REM - Typical  Overlay  Request 

30 1 60  F$= " ©AEROSOL /P R OGR AM / TASK 1 ” 

30170  GO  TO  30060 
30180  REM  i 

30190  REM 

30200  REM  - 

30210  REM - LOAD  OVL  AREA 

30220  IF  FS=Z*  THEN  30280 
30230  DELETE  20020, 29990 
30240  G=MEMORY 
30250  CALL  "UNIT",U 
30260  APPEND  F$; 20010, 10 
30270  2%=F% 

30280  RETURN 

20000  REM 

20010  REM  ! - 

20020  REM  I  © AER OSOL/ PROGRAM/ TASK 1 
20030  REM  !  STUART  GATHMAN  NRL  CODE  4327 

20040  REM  I - 

20050  REM 
20060  UNIT  1 
20070  DIM  A*<300) 

20080  PRINT  "IS  THIS  A  FRESH  START  <Y  OR  N>?"; 
20090  INPUT  B% 

20100  IF  "N"  THEN  20310 

20110  REM  SET  UP  NEW  PARAMETER  STORAGE  FILE 
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20120  GOSUB  22330 

20130  IF  B*="N"  THEN  20080 

20140  REM  SET  UP  NEW  OUTPUT  FILE 

20150  GOSUB  22490 

20160  IF  B*="N"  THEN  20080 

20170  REM  INITIATE  OUTPUT  FILE 

20180  GOSUB  22630 

20190  REM  DIMENSION  VARIABLES 

20200  GOSUB  20930 

20210  P=0 

20220  REM  INPUT  NON  MET  PARAMETERS 
20230  GOSUB  22750 

20240  REM  INPUT  Z  GRID  AND  MET  PARAMETERS  WITH  N1S<01  =1E“5 
20250  ES=MYM 
20260  GOSUB  22960 

20270  F$=" ©AEROSOL /PR QGR AM/ TASK2" 

20280  GO  TO  30060 
20290  REM 
20300  REM 

20310  PRINT  "IS  THIS  A  DUPLICATE  RUN  WITH  DIFFERENT  DRY  SIZE" 
20320  PRINT  ”<Y  OR  N  )  ?"; 

20330  INPUT  B* 

20340  IF  B*="N"  THEN  20560 

20350  REM  REOPEN  OLD  PARAMETER  STORAGE 

20360  GOSUB  22420 

20370  REM  OPEN  NEW  OUTPUT  FILE 

20380  GOSUB  22490 

20390  REM  INITIATE  OUTPUT  FILE 

20400  GOSUB  22630 

20410  IF  B*="N"  THEN  20080 

20420  REM  DIMENSION  VARIABLES 

20430  GOSUB  20930 

20440  READ  #2:  P,  N1 ,  K,  H9,  U9,  Z,  A 1 ,  C  1 ,  D 1 ,  SO,  RB,  V 

20450  REM  INPUT  NON  MET  PARAMETERS 

20460  GOSUB  22750 

20470  Nl  =  l.  OE-5 

20480  01=N1 

20485  P( 18>=0 

20490  P ( 16 ) =0 

20500  P ( 23 ) =1 

20510  P ( 24 ) =0 

20520  F*=,,©AER0S0L/PR0GRAM/TASK2" 

20530  GO  TO  30060 
20540  REM 
20550  REM 

20560  PRINT  "IS  THIS  A  SAVE  THE  WORLD  RECOVERY  (Y  OR  N)?"; 

20570  INPUT  B% 

20580  IF  B*="Y"  THEN  20820 

20590  PRINT  "AT  WHAT  FILE  NAME  IS  THE  INITIAL  CONCENTRATION  DATA’ 
20600  INPUT  B% 

20610  OPEN  B*; 3#  "R " #  A$ 

20620  PAGE 
20630  PRINT  A* 

20640  REM  OPEN  NEW  OUTPUT  DATA  FILE 
20650  GOSUB  22490 

20660  REM  INITIATE  NEW  OUTPUT  FILE 
20670  GOSUB  22630 

20680  REM  OPEN  NEW  PARAMETER  STORAGE  FILE 
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20690  GOSUB  22320 
20700  E**"N" 

20710  REM  DIMENSION  VARIABLES 
20720  QOSUB  20930 

20730  REM  INPUT  NON  MET  PARAMETERS 
20740  GOSUB  22750 

20750  REM  INPUT  1  GRID  AND  MET  MARAMETERS 
20760  GOSUB  22950 

20770  F$=  " ©AEROSOL /PROGR AM/ TASK2 11 
20780  GO  TO  30060 
20790  REM 
20800  REM 

20810  REM  RECOVER  FROM  A  SAVE  THE  WORLD 
20820  GOSUB  20930 

20830  OPEN  "©AER0S0L4/PARAMETER/ST0RAGE"; 2, «F"i AS 
20840  OPEN  ,,eAER0S0L4/0UTPUT“;  1,  "F",  A$ 

20850  READ  #2:  Pi  Nl,  K,  H9,  U9#  Z>  Al«  Cl#  Dl#  SO*  R8>  V 
20860  01=N1 

20870  F*=  ” ©AEROSOL /PROGRAM/ TASK 5 M 
20880  GO  TO  30060 
20890  REM 

20900  REM  ##############  SUBROUTINES  A####################### 

20910  REM 

20920  REM==========SUBROUTINE  D I  MENS I 0N= ==================== =========== 

20930  DIM  U9<24), H9<24), SO (24) 

20940  DIM  K ( 24 )  >  Z  <  24 ) *  01 <  24 ) *  R8 ( 24 ) *  V ( 24 ) 

20950  DIM  N1 (24), P<24), A1 <24)> D1 < 24 ) , Cl (24) 

20960  RETURN 
20970  REM 

20980  REM«====«====VERTICAL  SPACING  SUBROUT  I NE==="===============“= 

20990  REM 
21000  2  < 1 )  =  1 
21010  01  =  1.  OE-5 
21020  FOR  1=1  TO  23 
21030  Z < I +  1 > =Z < I ) +1 .  4^1+3.  5 
21040  01 ( 1  +  1 >«1.  0E-5 
21050  NEXT  I 
21060  Nl=01 

21070  IF  E$0"N"  THEN  21100 
21080  READ  #3,  144: N1 
21090  0 1 =N 1 
21100  RETURN 
21110  PRINT  Z 

21120  REM=====R0UT INE  TO  FILL  THE  INVERSION  HEIGHT  L I ST================ 

21130  REM 

21140  PRINT  ” INPUT  EXISTING  INVERSION  HEIGHT  (METERS)", 

21150  INPUT  H8 
21160  19=1 
21170  GOSUB  21270 

21180  PRINT  "ARE  THERE  FUTURE  ESTIMATES  OF  INVERSION  HEIGHTS  (V  OR  N)?M 
21190  INPUT  B% 

21200  IF  B*="Y"  THEN  21220 
21210  RETURN 

21220  PRINT  "INPUT  RELATIVE  HOUR  AND  INVERSION  HEIGHT" 

21230  INPUT  19,  H8 
21240  I 9= I NT ( 19 ) 

21250  GOSUB  21270 
21260  GO  TO  21180 
21270  FOR  1=19  TO  24 
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21280  H9< I  )=H8 
21290  NEXT  I 
21300  RETURN 
21310  REM 

21320  REM=====SUBROUTINE  TO  FILL  THE  WIND  SPEED  L I ST============== 

21330  REM 

21340  PRINT  "INPUT  EXISTING  10  METER  WIND  SPEED  <M/S>" 

21350  INPUT  H8 

21351  IF  H8>0  THEN  21360 

21352  HS=0.  1 
21360  19=1 
21370  GOSUB  21460 

21380  PRINT  "ARE  THERE  FUTURE  ESTIMATES  OF  THE  WIND  SPEED  (V  OR  N) 
21390  INPUT  B* 

21400  IF  B*="N"  THEN  21490 

21410  PRINT  "INPUT  NEXT  RELATIVE  HOUR  AND  WIND  SPEED" 

21420  INPUT  19,  HQ 

21421  IF  H8>0  THEN  21430 

21422  H8=0.  1 
21430  I 9= I NT ( 19 ) 

21440  GOSUB  21460 
21450  GO  TO  21380 
21460  FOR  1=19  TO  24 
21470  U9 ( I ) =H8 
21480  NEXT  I 

21490  RETURN 
21500  REM 

21510  REM==================  SST  INPUT  SUBROUT INE================== 

21520  REM 
21530  PRINT 

21540  PRINT  "INPUT  EXISTING  SEA  SURFACE  TEMPERATURE  VALUE  (CENT)" 

21550  INPUT  X 

21560  19=1 

21570  GOSUB  21660 

21580  PRINT  "ARE  THERE  FUTURE  ESTIMATES  OF  THIS  VALUE?  (V  OR  N)"i 
21590  INPUT  B* 

21600  IF  B*="N"  THEN  21650 

21610  PRINT  "INPUT  RELATIVE  HOUR  AND  DATA  VALUE"; 

21620  INPUT  19,  X 
21630  GOSUB  21660 
21640  GO  TO  21580 
21650  RETURN 
21660  FOR  1=19  TO  24 
21670  SO(I)=X 
21680  NEXT  I 
21690  RETURN 
21700  REM 

21710  REM==============AIR  TEMPERATURE  INPUT  SUBROUT I NE=========== 

21720  REM 
21730  PRINT 

21740  PRINT  "INPUT  EXISTING  AIR  TEMPERATURE  VALUE  (CENT)" 

21750  INPUT  X 
21760  19=1 
21770  GOSUB  21860 

21780  PRINT  "ARE  THERE  FUTURE  ESTIMATES  OF  THIS  DAT A^  (Y  OR  N>", 
21790  INPUT  B* 

21800  IF  B*="N"  THEN  21850 

21810  PRINT  "INPUT  RELATIVE  HOUR  AND  DATA  VALUE"; 

21820  INPUT  19,  X 
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21830  GOSUB  21860 
21840  GO  TO  21780 
21850  RETURN 
21860  FOR  1=19  TO  24 
21870  A1 ( I >=X 
21880  NEXT  I 
21890  RETURN 
21900  REN 

21910  REM=== ==========  DEWPOINT  INPUT  SUBROUT  I NE =====  =  ===  ===  =  =  =  === ^  -  === 

21920  REN 
21930  PRINT 

21940  PRINT  “INPUT  EXISTING  DEWPOINT  TENPERATURE  VALUE  (CENT)” 

21950  INPUT  X 
21960  19=1 
21970  GOSUB  22060 

21980  PRINT  “ARE  THERE  FUTURE  ESTINATES  OF  THIS  DATA?  (Y  OR  N) 

21990  INPUT  B* 

22000  IF  B*=MN“  THEN  22050 

22010  PRINT  “INPUT  RELATIVE  HOUR  AND  DATA  VALUE”, 

22020  INPUT  19, X 
22030  GOSUB  22060 
22040  GO  TO  21980 
22050  RETURN 
22060  FOR  1=19  TO  24 
22070  D 1 < I )=X 
22080  NEXT  I 
22090  RETURN 
22100  REN 

22110  REN  ==============  CLOUD  COVER  INPUT  SUBROUTINE  ========= 

22120  REN 
22130  PRINT 

22140  PRINT  “INPUT  EXISTING  CLOUDCO VER  IN  TENTHS*' 

22150  INPUT  X 
22160  19* 1 
22170  GOSUB  22260 

22180  PRINT  “ARE  THERE  FUTURE  ESTINATES  OF  CLOUD  COVER?" 

22190  INPUT  B* 

22200  IF  B*=“N“  THEN  22250 

22210  PRINT  “INPUT  RELATIVE  HOUR  AND  CLOUDCOVER  VALUE'* 

22220  INPUT  19,  X 
22230  GOSUB  22260 
22240  GO  TO  22180 
22250  RETURN 
22260  FOR  1*19  TO  24 
22270  Cl ( I )  =  X 
22280  NEXT  I 
22290  RETURN 
22300  REN 

22310  REN  ===SUBROUT INE  TO  SET  UP  NEW  PARAMETER  STORAGE  ====== 

22320  REN 

22330  CALL  “FILE",  1 ,  “@AER0S0L4/PARAMETER/ST0RAGE“ ,  A* 

22340  IF  LEN<AS)=0  THEN  22420 

22350  PRINT  “OK  TO  DESTROY  DATA  ON  ©AER0S0L4 /PARAMETER /STORAGE 
22360  PRINT  “TYPE  IN  'Y '  OR  'N'“, 

22370  INPUT  B* 

22380  IF  B*=“N“  THEN  22450 
22390  KILL  “ &AER0S0L4/P ARAMETER /STORAGE “ 

22400  CREATE  “©AER0S0L4/PARAMETER/ST0RAGE" ; 700,  0 
22410  REM  REOPEN  PARAMETER  STORAGE 
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22420  OPEN  ” ©AEROSOL4/PARAMETER/STOR AGE  '  .  2/  "F">  A* 

22430  PAGE 
22440  PRINT  A* 

22450  RETURN 
22460  REM 

22470  REM  ======SUBR0UT I NE  TO  SET  UP  NEW  OUTPUT  FILE  ======== 

22480  REM 

22490  CALL  "FILE " , 1 , "GAERQSOLA/OUTPUT" , AS 
22500  IF  LEN<A*)=0  THEN  22550 

22510  PRINT  "OK  TO  DESTROY  DATA  ON  ©AER0S0L4/ OUTPUT  (Y  OR  N)?“, 
22520  INPUT  B* 

22530  IF  BS= "N"  THEN  22590 
22540  KILL  "@AER0S0L4/0UTPUT" 

22550  CREATE  "©AER0S0L4/0UTPUT " .  "U" ;  240,  250 
22560  OPEN  "@AER0S0L4/0UTPUT“ ; 1, "F"» A« 

22570  PAGE 
22580  PRINT  A* 

22590  RETURN 
22600  REM 

22610  REM  ==  SUBROUTINE  TO  INITIATE  RANDOM  DISC  FILE  #1  ====== 

22620  REM 
22630  DIM  DS(657) 

22640  D*=" " 

22650  FOR  1=1  TO  245 
22660  DS=D$$<"  " 

22670  NEXT  I 

22680  FOR  1=1  TO  240 

22690  WRITE  #1,  I  :  D% 

22700  NEXT  I 
22710  RETURN 
22720  REM 

22730  REM  =====SUBROUT I NE  TO  INPUT  NON  MET  P AR AMETERS========= 

22740  REM 
22750  P  <  22 )  =  1 
22760  P  (  1 9  )  =0.  1 
22770  P ( 20 ) = 1 

22780  PRINT  “AT  WHAT  LATITUDE  IS  THIS  TAKING  PLACE?” 

22790  INPUT  P<5) 

22800  PRINT  "ENTER  THE  DIAMETER  OF  THE  DROPLETS  IN  MICRONS" 
22810  INPUT  P( 1 ) 

22820  PRINT  "ENTER  THE  REFERENCE  HUMIDITY  OF  DROPLET  SIZE" 

22830  INPUT  P(2) 

22840  IF  P <2X99.  9  THEN  22860 
22850  P  (  2  )  =99.  9 

22860  P<3)=P  <  1  )*<  1+0.  9/  <  1  -P  <  2 )  / 1 00 )  )**'<- 1  /  3 ) 

22870  IF  P ( 3 ) >90  THEN  22900 
22880  IF  P  (3X0  07  THEN  22900 
22890  RETURN 

22900  PRINT  "DROPLET  SIZE  IS  OUTSIDE  OF  MODEL  LIMIT:  TRY  AGAIN" 
22910  GO  TO  22800 
22920  REM 

22930  REM  =====  SUBROUTINE  TO  INPUT  Z  GRID  AND  MET  PARAMETERS  = 
22940  REM 

22950  REM  WORK  OUT  VERTICAL  SPACING 
22960  GOSUB  21000 
22970  H9=0 

22980  REM  INPUT  INVERSION  HEIGHT 
22990  GOSUB  21140 
23000  REM 
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23010  REM  INPUT  WIND  SPEED 
23020  GOSUB  21340 
23030  REM 

23040  REM  INPUT  SEA  SURFACE  TEMPERATURE  VALUES 
23050  GOSUB  21530 
23060  REM 

23070  REM  INPUT  AIR  TEMPERATURE  VALUES 
23080  GOSUB  21730 
23090  REM 

23100  REM  INPUT  DEWPOINT  TEMPERATURE  VALUES 
23110  GOSUB  21930 
23120  REM 

23130  REM  INPUT  CLOUD  COVER  DATA 
23140  GOSUB  22130 
23150  RETURN 
20000  REM 

20010  REM  ; - I 

20020  REM  ,i  ©AEROSOL / PROGRAM /TASK2  S 

20030  REM  {  S  GATHMAN  NRL  CODE  4327  ! 

20040  REM  { - ; 

20050  REM 
20060  REM 

20070  REM==ROUT INE  TO  UPDATE  THE  METEOROLOGICAL  TERMS  OF  P  VECTOR== 
20080  REM  CALCULATE  DRAG  COEFFICIENT#  CIO  (ROLL  FIG  41  P160) 

20090  P ( 6 ) =0.  0011+4  0E-5*U9 ( P ( 22 )  ) 

20100  REM  CALCULATE  THE  FRICTION  VELOCITY, U*  (M/SEC) (ROLL  EQN  4  51) 
20110  P ( 8 ) =U9 ( P ( 22 ) ) *SQR ( P  <  6 )  ) 

20120  REM  CALCULATE  ZZERO(M)  CROLL  EQN  4  183 
20130  P(9)=0.  00 1 25*P ( 8 ) ' 2 
20140  REM  ROLL  EQN  4.  14 

20150  REM  CALCULATE  BULK  RICHARDSON'S  NUMBER, RBI 
20160  T=S0 ( P ( 22 ) ) 

20170  GOSUB  21260 
20180  VO=E 
20190  T =D 1 (P (22) ) 

20200  GOSUB  21260 
20210  V 1 =E 

20220  R0=V0*0.  622*0  98/1013  25 
20230  R=V1*0.  622/1013  25 

20240  T=(A1 (P(22) )+273)*( l+R/O. 622) /(1+R) 

20250  T0= ( SO ( P ( 22 ) )+273)*< l+RO/O.  622)/<l+R0) 

20260  T =T * 1 .  0098 

20270  P( 10)=(T-T0)*3.  55/U9< P ( 22 ) ) A2 
20280  REM 

20290  REM  CALCULATE  THE  Z/L  PARAMETER  AT  10  METERS  B ARKER&B AXTER ( 75 ) 
20300  REM 

20310  P<4>*LOG< 10/P<9> I/O.  4 
20320  IF  P  <  4 ) => 1 0  THEN  20390 
20330  RESTORE  20340 

20340  DATA  0.  07872,  0.  05532,  0.  006197,  4  7,  0  35 
20350  READ  CO,  C2,  C3,  C4,  C5 

20360  P(7)*C5*P(4)*(P( 1 0 ) -CO+SQR ( C2*P ( 10)+C3) )/( 1+C4*P( 10) ) 

20370  P ( 7 ) -P  <  7 ) / 1 0 
20380  GO  TO  20450 

20390  P  <  7 ) =P ( 10 ) * ( 0.  47 1 *P  <  4 ) - 1 .  045) 

20400  IF  P ( 7 ) >— 0.  05  THEN  20330 
20410  REM 

20420  REM  UPDATE  THE  FRICTION  MIXING  RATIO, Q* 
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20430  REM  REQUIRES  P ( 6 ) , P < 22 ) , SO, D 1 
20440  REM  RESETS  P<12) 

20450  Zl  =  10 
20460  T =D 1  (P (22)  ) 

20470  GOSUB  21360 
20480  Q 1 =R  5 
20490  Z1=0 
20500  T=SO(P (22) ) 

20510  GOSUB  21360 
20520  G0=0  98*R5 

20530  P( 1 2 ) =SGR ( P  <  6 ) )*<Ql-Q0)/0  38 
20540  REM 

20550  REM  UPDATE  THE  FRICTION  POTENTIAL  TEMPERATURE, THETA* 
20560  REM 

20570  REM  REQUIRES  P < 6 ) ,  P ( 22 ) ,  SO,  A  1 
20580  REM  RESETS  P ( 1 1 ) 

20590  Zl  =  10 
20600  T -A 1 (P (22) ) 

20610  GOSUB  21430 
20620  T 1 =T9 
20630  Z1=0 
20640  T =S0 ( P ( 22 ) ) 

20650  GOSUB  21430 

20660  P( 1 1 ) =SQR ( P  <  6 ) )*<Tl-T0)/0  38 
20670  REM 

20680  REM  UPDATE  THE  HEIGHT  OF  SUPERADI OB AT I C  ZONE, ZSA 
20690  REM  REQUIRES  P(22),S0, A1 
20700  REM  RESETS  P(13) 

20710  D=ABS  <  A1 ( P ( 22 )  ) -SO ( P ( 22 )  ) ) 

20720  IF  D>0.  5  THEN  20750 

20730  P( 13)=10*D 

20740  GO  TO  20810 

20750  P( 13)  =  11  92+9.  69*L0G<D) 

20760  REM 

20770  REM  UPDATE  THE  HEIGHT  OF  THE  CONDENSATION  LAYER 
20780  REM  REQUIRES  ZSA ( P ( 1 3 ) ) ,  AND  POTENTIAL  TEMP 
20790  REM  AND  MIXING  RATIO  AT  TOP  OF  SUPER AD I OBAT IC  ZONE 
20800  REM  RESETS  P(14) 

20810  Z1=P( 13) 

20820  GOSUB  21650 
20830  Z 1=P ( 13) 

20840  GOSUB  21760 
20850  T=T7 
20860  R=Q7 
20870  Z 1=P  < 13 ) 

20880  GOSUB  21500 
20885  T6=T6+273.  15 
20890  GOSUB  21580 
20900  X=EXP < LOG <T6)/0.  286) 

20910  XI =732.  02-150.  4 1  * < LGT < X ) -LGT < V6 ) ) 

20920  X 1  =  X 1 +7.  2 1  *  <  LGT ( X ) -LGT  <  V6 ) >^2 
20930  X1  =  X  1+273. 15 
20935  T7=T7+273.  15 

20940  P 1  =  1000*EXP  < -LOG  <  T7/X 1 )/0.  286) 

20950  GOSUB  21050 

20960  P( 14)=A 

20961  IF  A>0 .  3*H9 ( P ( 22 ) )  THEN  20970 

20962  P( 14)=H9<P(22) ) 

20970  F*=" ©AEROSOL /PROGR AM/ TASK3" 
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20980  GO  TO  30060 

20990  REtl##################  SUBROUTINES  ######################### 

21000  REN 
21010  REN 

21020  REN  ****************  ALT  P  1  )  -  A  **************************** 

21030  REN 

21040  REN  INPUT  P ( MB )  OUTPUT  A  <  METERS ) 

21050  IF  P 1 <» 1013  THEN  21080 

21060  A-0 
21070  RETURN 

21080  IF  PI <=958  THEN  21110 
21090  A-9  09* (101 3-P 1 ) 

21100  RETURN 

21110  A=7Q50*L0G( 1021  38/P1) 

21120  RETURN 
21130  REM 

21140  REM  ***************  PRESSURE  £(Z=A)  *************************** 

21150  REN 

21X60  REM  INPUT  A ( METERS)  OUTPUT  PI (MB) 

21170  IF  A <500  THEN  21200 

21180  P 1 = 1 02 1  38*EXP ( - 1  2739E-4*A) 

21190  RETURN 

21200  P 1  =  1 01 3— 55*A/  500 

21210  RETURN 

21220  REM 

21230  REM  **************  £=VAPPR(T>  ************************* 

21240  REM 

21250  REM  INPUT  T ( DEGREES  C)  OUTPUT  E ( MB ) 

21260  Tl  =  l~373.  25/ (T+273  25) 

21270  RESTORE  21280 

21280  DATA  1013.  25,  13.  3185,  1  976,  0.  6445,  0  1299 
21290  READ  RO,  R  1 ,  R2,  R3,  R4 

21300  E=RO*EXP  (R1*T  1 -R2*T 1  'S2'-R 3*T1''3-R4*T1"'4) 

21310  RETURN 
21320  REM 

21330  REM  **********  R 5=SMI XR ( Z 1 , T )  ************************* 

21340  REM 

21350  REM  INPUT  Z1  (METERS)  , T ( DEGREES  C)  OUTPUT  R5  (GN/KG) 

21360  GOSUB  21260 
21365  A=Z 1 
21370  GOSUB  21170 
21380  R5=0.  622*E/(P1-E) 

21390  RETURN 
21400  REM 

21410  REM  ***  CONVERTS  TEMP  TO  POTENTIAL  TEMP  G  ALT  Z1  ********** 

21420  REM  INPUT  T ( DEGREES  C ),  Z 1 ( METERS )  OUTPUT  POTENTIAL  TEMP  T ( DEG  C) 
21430  GOSUB  21170 

21440  T9~  <  T+273.  1 5 > *EXP < O.  286*L0G ( 1 000/P  1 > > 

21450  RETURN 
21460  REM 

21470  REM  ***  CONVERTS  POTENTIAL  TEMP  G  ALT  Z1  TO  TEMPERATURE  **** 

21480  REM 

21490  REM  INPUT  T  <  DEG  C)  $<  Z  1  <  METERS ) :  OUTPUT  T6  (  DEG  C) 

21500  GOSUB  21170 

21510  T6=(T+273.  15)/EXP(0.  286*L0G (  1000/P  1  )  ) 

21520  T6=T6-273.  15 
21530  RETURN 
21540  REM 

21550  REM  ******CONVERTS  MIXING  RATIO  G  ALT  Z1  TO  VAPOR  PRESSURE 
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21560  REM 

21570  REM  INPUT  R  <  G/KG  )  ZtZ  1  ( METERS  )  OUTPUT  V6 ( MB ) 

21580  GOSUB  21170 
21590  V6=Pl*R/<0  622+R) 

21600  RETURN 
21610  REM 

21620  REM  ******  CALCULATION  OF  MIXING  RATIO  IN  SUPER AD I OBAT I C  ZONE  ** 
21630  REM 

21640  REM  INPUT  Z 1 , SO, P < 22 ) > P < 1 2 ) , P < 9 ) , P ( 7 )  OUTPUTS  07 
21650  X=<Z1+P<9) )/P<9) 

21660  Z1=0 
21670  T=S0(P (22) ) 

21680  GOSUB  21360 
21690  G0=0  98*R 5 

21700  Q7=Q0+P  (  1  2  )  *  (  LOG  (  X  )  +4 .  8*P(7)  ) 

21710  RETURN 
21720  REM 

21730  REM  SUBROUTINE  TO  CALCULATE  POT  TEMP  IN  SUPERAD IOBAT I C  ZONE 
21740  REM 

21750  REM  INPUT  Z 1 , A 1 , P ( 9 ) , P ( 22 ) , P ( 1 1 ) > P ( 7 ) :  OUTPUT  17 
21760  X=<  Zl+P<9) >/P<9> 

21770  Z1=0 
21780  T=A1 (P (22) ) 

21790  GOSUB  21500 

21800  T7=T6+P ( 1 1 )*<LOG<  X ) +4  8*P (7 ) ) 

21810  RETURN 

20000  REM 

20010  REM  ! - ! 

20020  REM  !  ©AEROSOL /PROGR AM /TASK3  ! 

20030  REM  I  AFTER  NRL  REPORT  8279  :  S  GATHMAN  \ 

20040  REM  ! - ; 

20050  REM 

20060  REM  ROUTINE  TO  CALCULATE  REL  HUM  PROFILE  FROM  P  VALUES 
20070  REM 

20080  IF  P  <  7 ) >0.  1  THEN  20120 
20090  IF  P ( 7 ) <— 0.  1  THEN  20140 
20100  IF  Cl (P(22> )>2.  5  THEN  21390 
20110  GO  TO  21690 
20120  IF  C 1 ( P ( 22 ) ) >1  THEN  20200 
20130  00  TO  20790 

20140  IF  Cl <P(22) )>2.  5  THEN  22370 
20150  GO  TO  22630 
20160  REM 

20170  REM - 

20180  REM 

20190  REM  CASE  1  (STABLE  WITH  C/C  >  107.) 

20200  FOR  1=1  TO  24 

20210  IF  Z ( I ) >20  THEN  20490 

20220  REM  FOR  ALTITUDES  BELOW  20  METERS 

20230  Z1=0 

20240  T=SO(P (22) ) 

20250  GOSUB  23600 
20260  T 0=T9 
20270  Z 1=20 
20280  T=A1 ( P ( 22 ) ) 

20290  GOSUB  23600 
20300  T1=T9 

20310  T7=T0+Z< I )*<Tl-T0)/20 
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20320  T=D1 ( P  <  22 ) ) 

20330  GOSUB  23420 
20340  A=20 
20350  GOSUB  23330 
20360  R6=622*E/ (PI -E ) 

20380  T=S0 ( P ( 22 ) ) 

20390  GOSUB  23420 
20400  E=0.  98*E 
20410  A=0 
20420  GOSUB  23330 
20430  R0=622*E/ (Pl-E) 

20440  Q7=R0+Z( I )*<R6-R0>/20 
20450  Z1=Z(I) 

20470  GOSUB  24100 
20480  GO  TO  20730 

20490  IF  Z  < I > >H9 ( P ( 22 ) )  THEN  20720 
20500  REM  BELOW  STRATUS  CLOUD 
20510  Z 1=20 
20520  T=A1 <  P ( 22 ) ) 

20530  GOSUB  23600 

20540  T7=T9+2* ( Z ( I > -20 ) / ( H9 ( P < 22 ) ) -20) 

20550  T =D 1 ( P ( 22 ) ) 

20560  Z 1=20 
20570  GOSUB  23520 
20590  Q0=R5*1000 
20600  T=T9+2 
20610  Z 1 =H9 ( P ( 22  > ) 

20620  GOSUB  23690 
20630  T=T6 
20640  Z 1=H9  <  P  <  22 ) ) 

20650  GOSUB  23520 
20660  Q1=R5*1000 

20670  Q7=Q0+(Z( I ) -20) * ( 01 -QO ) / ( H9 < P ( 22 ) )-20) 
20680  Z1=Z(I> 

20700  GOSUB  24100 
20710  GO  TO  20730 
20720  R8  < I ) =100 
20730  NEXT  I 
20740  GO  TO  23140 
20750  REM 

20760  REM - 

20770  REM 

20780  REM  CASE  2  (STABLE  WITH  C/C  <  107.) 

20790  FOR  1=1  TO  24 

20800  IF  Z ( I ) >20  THEN  21060 

20810  REM  FOR  ALTITUDES  BELOW  20  M 

20820  Z 1=0 

20830  T=SO  <  P ( 22 ) ) 

20840  GOSUB  23600 
20850  T0=T9 
20860  Z 1 =20 
20870  T=A1 <  P  <  22  > ) 

20880  GOSUB  23600 
20890  T1=T9 

20900  T7=T0+Z ( I ) * ( Tl-TO ) /20 
20910  T=D1 ( P ( 22 ) ) 

20920  GOSUB  23420 
20930  A*20 
20940  GOSUB  23330 
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20950  R6=622*E/<P1-E) 

20960  T=S0 ( P  <  22 ) ) 

20970  GOSUB  23420 
20980  E=0.  98*E 
20990  A=0 
21000  GOSUB  23330 
21010  R0=622*E/(P1-E) 

21020  Q7=R0+Z(I)*(R6-R0)/20 
21030  Z1=Z(I> 

21040  GOSUB  24100 
21050  GO  TO  21330 

21060  IF  Z  < I ) >H9 ( P ( 22 ) )  THEN  21190 

21070  REh  FOR  ALTITUDES  BELOW  INVERSION  LEVEL 

21080  Z 1 =20 

21090  T=A1 ( P l 22 ) ) 

21100  GOSUB  23600 

21110  T7=T9+2*  <  Z  <  I  )  -20 )  /  ( H9  <  P  ( 22 )  )  -20 ) 

21120  T=D1 ( P ( 22 ) ) 

21130  Z  15=20 
21140  GOSUB  23520 
21150  G7=R 5*1000 
21160  Z1=Z(I) 

21170  GOSUB  24100 
21180  GO  TO  21330 

21190  IF  Z< I )>H9(P<22> )+500  THEN  21320 

21200  REM  ALTITUDES  IN  500  M  LAYER  ABOVE  INVERSION -BASE 
21210  Z 1=20 
21220  T=A 1 ( P ( 22 ) ) 

21230  GOSUB  23600 

21240  T7=T9+2+3*  <  Z ( I ) -H9 ( P ( 22 ) ) )/500 
21250  T=D1 (P (22) ) 

21260  Z 1=20 
21270  COSUB  23520 
21280  Q7=R5*1000 
21290  Z1=Z<I) 

21300  GOSUB  24100 
21310  GO  TO  21330 
21320  R8( I ) =78 
21330  NEXT  I 
21340  GO  TO  23140 
21350  REM 

21360  REM - 

21370  REM 

21380  REM  CASE  3  (NEUTRAL  WITH  C/C  >  257.) 

21390  FOR  1=1  TO  24 

21400  IF  Z ( I ) >20  THEN  21500 

21410  REM  FOR  ALTITUDES  BELOW  20  METERS 

21420  T=D1 <P ( 22 ) ) 

21430  GOSUB  23420 
21440  E1=E 
21450  T«A1 (P ( 22 ) ) 

21460  GOSUB  23420 

21470  R1=100*E1/E 

21480  RB  < I ) =98+ ( R 1 -98 ) *Z  < I ) /20 

21490  GO  TO  21660 

21500  IF  Z  < I ) >P  < 14 )  THEN  21650 

21510  REM  FOR  ALTITUDES  BELOW  ZCON 

21520  T=D1 <P (22) ) 

21530  GOSUB  23420 
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21540  Z 1 =20 
21550  COSUB  23330 
21560  07*0.  622*E/<P1HE) 

21570  T=SO(P (22) ) 

21500  Z1=0 
21590  GOSUB  23600 
21600  T7=S0 ( P  <  22 ) ) 

21610  Z1=Z(I> 

21620  07=1000*07 
21630  GOSUB  24100 
21640  GO  TO  21660 
21650  R8< I ) =100 
21660  NEXT  1 
21670  GO  TO  23140 

21680  REM  CASE  4  (NEUTRAL  WITH  C/C<25%> 

21690  FOR  1=1  TO  24 

21700  IF  Z ( I ) >20  THEN  21800 

21710  REM  FOR  ALTITUDES  BELOW  20  METERS 

21720  T=D1 <  P ( 22 ) ) 

21730  GOSUB  23420 
21740  El -E 
21750  T=A1 (P (22) ) 

21760  GOSUB  23420 
21770  R1=100*E1/E 
21780  R8(I)=98+(Rl-98)*Z(I>/20 
21790  GO  TO  22310 

21800  IF  Z(I)>0.  8*PU4)  THEN  21960 

21810  REM  FOR  ALTITUDES  IN  HOMOGENEOUS  LAYER 

21820  T=D1 (P (22) ) 

21830  GOSUB  23420 
21840  Z 1=20 
21850  GOSUB  23330 
21860  G7=0.  622*E/(P1-E) 

21070  T=S0(P (22) ) 

21880  Z1=0 
21890  GOSUB  23600 
21900  T7=T9 
21910  T7=S0(P (22) ) 

21920  Z1*Z(I> 

21930  07=07*1000 

21940  GOSUB  24090 

21950  GO  TO  22310 

21960  IF  Z(I)>P(14)  THEN  22300 

21970  REM  FOR  ALTITUDES  IN  THE  TRANSITION  ZONE 

21980  Z5*Z(I)-P(14)*0  8 

21990  Z1=0 

22000  T=SO(P (22) ) 

22010  GOSUB  23600 
22020  Z1=0.  8*P(14) 

22030  T*T9 
22040  GOSUB  23690 
22050  T»T6+0.  006*Z5 
22060  Z1*Z(I> 

22070  GOSUB  23600 
22080  T7*T9 

22090  REM  FIND  MIXING  RATIO  AT  ZCON 
22100  Z 1 =P ( 14) 

22110  T=T6+0.  006*0.  2*P<14) 

22120  GOSUB  23600 
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22130  T=T9 
22140  GOSUB  23420 
22150  VI =0.  78*E 
22160  Z 1  =P  ( 14) 

22170  GOSUB  23330 

22180  M9=0  622*V1/(P1-V1 ) 

22190  REM  FIND  MIXING  RATIO  AT  8*ZC0N 
22200  T =D1 (PC22) ) 

22210  GOSUB  23420 
22220  21=20 
22230  GOSUB  23330 
22240  M8=0.  622*E/(P1-E> 

22250  X=M8+ < M9-M8 ) * C  2 ( I ) -0.  8*P(14))/<0.  2*P(14)> 

22260  Q7=X*1000 
22270  2 1  =  Z  < I ) 

22280  GOSUB  24100 
22290  GO  TO  22310 
22300  R8 ( I ) =78 
22310  NEXT  I 
22320  GO  TO  23140 
22330  REM 

22340  REM - 

22350  REM 

22360  REM  CASE  5  (UNSTABLE  WITH  C/C  >  257.) 

22370  FOR  1=1  TO  24 
22380  IF  Z ( I ) >P  < 13  >  THEN  22470 
22390  REM  FOR  SUPERADIOBADIC  LAYER 
22400  Z1=Z<I) 

22410  GOSUB  23850 
22420  Z1=Z(I) 

22430  GOSUB  23960 
22440  Z1=Z<I) 

22450  GOSUB  24100 

22460  GO  TO  22570 

22470  IF  Z ( I ) >P (14)  THEN  22560 

22480  REM  FOR  ALTITUDES  BETWEEN  ZSA  AND  ZCON 

22490  Z 1 =P \ 13) 

22500  GOSUB  23850 
22510  Z 1=P ( 13 ) 

22520  GOSUB  23960 
22530  Z1=Z(I) 

22540  GOSUB  24100 
22550  GO  TO  22570 
22560  R8  <  I )  =  1 00 
22570  NEXT  I 
22580  GO  TO  23140 
22590  REM 

22600  REM - 

22610  REM 

22620  REM  **  CASE  6  (UNSTABLE  ATMOSPHERE  WITH  C/C  <  257.) 

22630  FOR  1=1  TO  24 

22640  IF  Z  < I ) >P (13)  THEN  22730 

22650  REM  FOR  ALTITUDES  <  ZSA 

22660  Z1«Z(I) 

22670  GOSUB  23850 
22680  Z1*Z(I> 

22690  GOSUB  23960 
22700  Z1»Z(I) 

22710  GOSUB  24100 
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22720  GO  TO  23130 

22730  IF  Z(I>>0.  8*P<14>  THEN  22820 

22740  REM  FOR  ALTITUDES  IN  HOMOGENEOUS  ZONE 

22750  Z  1=P  < 13) 

22760  GOSUB  23850 
22770  Z  1=P  ( 13) 

22780  GOSUB  23960 
22790  Z1=Z < I ) 

22800  GOSUB  24100 

22810  GO  TO  23130 

22820  IF  Z  < I ) >P (14)  THEN  23120 

22830  REM  FOR  ALTITUDES  IN  TRANSITION  ZONE 

22840  Z5=Z< I )-P< 14)*0.  8 

22850  Z 1 =P ( 13 ) 

22860  GOSUB  23960 
22870  Z1=P< 14 ) *0.  8 
22880  T=T7 
22890  GOSUB  23690 
22900  T=T6+0.  006*Z5 
22910  Z1=Z(I) 

22920  GOSUB  23600 
22930  T7=T9 

22940  REM  FIND  MIXING  RATIO  AT  Z(I) 

22950  Z 1=P (14) 

22960  T-T6+0.  006*0.  2*P  < 14) 

22970  GOSUB  23600 
22980  T =T9 
22990  GOSUB  23420 
23000  V1=0.  78*E 
23010  A=P ( 14 ) 

23020  GOSUB  23330 
23030  R9=62?*V1 / <  P 1 -VI ) 

23040  Z 1»P (13) 

23050  GOSUB  23850 

23060  X3SG7+  ( R9-G7 )  *  <  Z  <  I )  -0.  8*P<14)  )/<0.  2*P<14>  ) 

23070  Q7=X 
23080  Z1=Z<I> 

23090  GOSUB  24100 

23100  GO  TO  23130 

23110  REM  CLOUDLESS  CLOUD  LAVER 

23120  R8 ( I ) =78 

23130  NEXT  I 

23140  F$=n ® AEROSOL /PR0GRAM/TASK4" 

23150  GO  TO  30060 

23160  REM  ##################  SUBROUTINES  ################*#########*###* 
23170  REM 

23180  REM  ****************  ALT(P1)-A  **************************** 

23190  REM 

23200  REM  INPUT  P(MB)  :  OUTPUT  A (METERS) 

23210  IF  Pl<=1013  THEN  23240 
23220  A*0 
23230  RETURN 

23240  IF  PI <*958  THEN  23270 
23250  A-9  09*< 1013-PI) 

23260  RETURN 

23270  A*7850*L0G< 1021  38/P 1) 

23280  RETURN 
23290  REM 

23300  REM  ***************  PRESSURE  @<Z=A)  *************************** 
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23310  REM 

23320  REM  INPUT  A ( METERS )  :  OUTPUT  P1<MB> 

23330  IF  AC500  THEN  23360 

23340  P I = 1 02 1  38*EXP<~1  2739E-4*A> 

23350  RETURN 
23360  Pl=1013~55*A/50 0 
23370  RETURN 
23380  REM 

23390  REM  **************  E=VAPPR(T)  *******  *-r-**********«-*-&r-** 

23400  REM 

23410  REM  INPUT  T  <  DEGREES  C)  OUTPUT  E ( MB ) 

23420  T  1  =  1-373.  25/ <  T+273.  1  5  ) 

23430  RESTORE  23440 

23440  DATA  1013.  25.  13.  3185,  1.  976,  0  6445,  0.  1299 
23450  READ  RO,  R  1 ,  R2,  R3,  R4 

23460  E^RO*EXP  <R1*T1-R2*T1  ''2-R3*Tl  '3-R4*Tl  "4  ) 

23470  RETURN 
23480  REM 

23490  REM  **********  R 5=SM I  X R < Z 1 ,  T )  ************************* 

23500  REM 

23510  REM  INPUT  Z1  (METERS)  , T ( DEGREES  C)  : OUTPUT  R5  (GM/GM) 

23520  GOSUB  23420 
23530  A=Z 1 
23540  GOSUB  23330 
23550  R5=0.  622*E/ ( P 1 -E ) 

23560  RETURN 
23570  REM 

23580  REM  ***  CONVERTS  TEMP  TO  POTENTIAL  TEMP  ©  ALT  Z1  ********** 
23590  REM  INPUT  T ( DEGREES  C ), Z 1 ( METERS )  : OUTPUT  POTENTIAL  TEMP  T ( DEG  C) 

23600  A=Z 1 
23610  GOSUB  23330 

23620  T9=  <  T+273.  15>*EXP (0.  286*L0G< 1000/PI )  ) 

23630  T9=T9~273. 15 
23640  RETURN 
23650  REM 

23660  REM  ***  CONVERTS  POTENTIAL  TEMP  @  ALT  Z1  TO  TEMPERATURE  **** 
23670  REM 

23680  REM  INPUT  T ( DEG  C)  &  Z 1 < METERS ) :  OUTPUT  T6 ( DEG  C) 

23690  A=Z 1 
23700  GOSUB  23330 

23710  T6=<T+273.  1 5 ) /EXP < O.  286*L0G < 1000/P  1 ) ) 

23720  T6=T6-273  15 
23730  RETURN 
23740  REM 

23750  REM  ******C0NVERTS  MIXING  RATIO  ®  ALT  Z1  TO  VAPOR  PRESSURE 
23760  REM 

23770  REM  INPUT  R  <  G/KG )  &Z1  (METERS)  :  OUTPUT  V6(MB) 

23780  GOSUB  23330 
23790  V6=Pl*R/<622+R > 

23800  RETURN 
23810  REM 

23820  REM  SUBROUTINE  TO  CALCULATE  MIXING  RATIO  IN  SUPERADIOBAT I C  ZONE 
23830  REM 

23840  REM  INPUT  Z 1 , SO*  P ( 22 )»P<12)*P(9)#P(7>:  OUTPUTS  G7 
23850  X*  <  Z 1 +P  <  9 ) )/P(9) 

23860  Z1=0 
23870  T»S0<  P  <  22 ) ) 

23880  GOSUB  23520 
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23890 

23900 

23910 

23920 

23930 

23940 

23950 

23960 

23970 

23980 

23990 

24000 

24010 

24020 

24030 

24040 

24050 

24060 

24070 

24080 

24090 

24100 

24110 

24120 

24130 

24140 

24150 

24160 

24170 

24180 

24190 

20000 

20010 

20020 

20030 

20040 

20050 

20060 

20070 

20071 

2007 2 

20073 

20074 
20080 
20090 
20100 
201  10 
20120 
20130 
20140 
20150 
20160 
20170 
20180 
20190 
20200 
20210 
20220 


Q0=980*R5 

Q7=Q0+P ( 1 2 ) *  1000* ( LOG ( X ) +4.  8*P (7) ) 

RETURN 

REM 

REM  SUBROUTINE  TO  CALCULATE  POT.  TEMP  IN  SUPER AD  I OB AT I C  ZONE 
REM 

REM  INPUT  Zl, Al, P (9), P (22). P( 1 1 )>  P(7)  :  OUTPUT  77 
X=<Z1+P<9> )/P<9) 

Z1=0 

T -A 1 ( P ( 22 ) > 

GOSUB  23690 

T7=T6+P  < 1 1 ) * ( LOG ( X ) +4  8*P  <  7 ) ) 

RETURN 

REM 

REM  SUBROUTINE  TO  CALCULATE  RELATIVE  HUMIDITY  FROM 

REN  INPUTS  OF  POTENTIAL  TEMPERATURE  (DEG  C),  ALT  (METERS) 

REM  AND  MIXING  RATIO  (G/KG)  AND  TO  OUTPUT  THIS 
REM  INFORMATION  IN  PERCENT  IN  VARIABLE  R8(I> 

REM 

REM  INPUTS  I , T7 *  07,  Z1 
REM  OUTPUTS  R8< I ) 

T ~T7 

GOSUB  23690 
R-Q7 

GOSUB  23780 
T=T6 

GOSUB  23420 

R8< I >=100*V6/E 

IF  R8(I>0100  THEN  24190 

R8 ( I ) - 1 00 

RETURN 

REM 

REM  ! - ! 

REM  J  ©AEROSOL /PROGRAM/ TASK4  ! 

REM  J  S.  GATHMAN  NRL  CODE  4327  ! 

REM  \ - I 

REM 

REM  CALCULATE  K  PROFILE  FROM  SURFACE  METEOROLOGY  DATA 
REM  K1=K<  1  )  S<  K3=K  '  (  1  ) 

P7=0.  1  *P ( 7  > 

IF  P(  10X1/4.  7  THEN  20080 
K=0 

GO  TO  20450 

IF  P(7K*0  THEN  20180 

REM 

REM - STABLE  CASE - 

REM 

K1  “O.  35*P<8)/<0  74+4.  7*P7) 

K3=0.  43*P<8>*P7/<0.  74+4.  7*P7X2 

GO  TO  20210 

REM 

REM - UNSTABLE  CASE - 

REM 

K 1 =0  47*P<8)*SGR< 1-9*P7) 

K3=Kl-2.  13*P(B)*P7/SQR< 1-9*P7) 

REM  CALCULATE  THE  K  PROFILE 
H8=H9 ( P ( 22 ) ) 

IF  P7O0  THEN  20250 


I 
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20230 

20240 

20250 

20260 

20270 

20280 

20290 

20300 

20310 

20320 

20330 

20340 

20350 

20360 

20370 

20380 

20390 

20400 

20410 

20420 

20430 

20440 

20450 

20460 

20470 

20480 

20490 

20500 

20510 

20520 

20530 

20540 

20550 

20560 

20570 

20580 

20590 

20610 

20620 

20630 

20640 

20650 

20660 

20680 

20690 

20700 

20705 

20710 

20720 

20730 

20740 

20750 

20760 

20770 

20780 

20790 

20800 

20810 


SET  DEGREES 

H8-0.  23*P<8)/<2*7  2921 1E-5*SIN< P < 5) ) ) 

K9=0 

F=K3+2*K1/(H8-Z< 1 ) ) 

F=F/ (H8-Z  < 1 ) ) ^2 
G=K 1 / ( HB-Z ( 1 ) )^2 
D=0 
K=0 

A=G-F*<2*H8+Z< 1 ) ) 

B=F*H8^2-2*H8* ( G-F*Z  < 1 ) ) 

FOR  1=1  TO  24 
IF  Z(IK=H8  THEN  20370 
K ( I ) =0 
GO  TO  20380 

K  < I >=F*Z< I ) "3+A*Z< I ) '2+B*Z< I ) +H8"2* ( G-F*Z ( 1 ) ) 

NEXT  I 

REM - 

REM  CALCULATE  THE  FALL  VELOCITY  PROFILE  VECTOR 
REM 

REM  INPUT  RB>  THE  RELATIVE  HUMIDITY  PROFILE  VECTOR 

REM  INPUT  P  <  3  >  >  THE  DRY  DROPLET  DIAMETER 

REM  OUTPUT  V9/  THE  FALL  VELOCITY  PROFILE  VECTOR 

FOR  1=1  TO  24 

IF  R8 ( I ) <99.  9  THEN  20480 

R8( I >=99  9 

D9=P<3)*< 1+0  9/ ( 1~R8< I >/100) >^( 1/3) 

V  < I ) =3  OE-5*  < 1+0  1 66/D9 ) *D9  '2 

NEXT  I 

REM 

REM  ROUTINE  TO  CALCULATE  THE  FLUX  OF  AEROSOL  FROM  THE  SURFACE  INTO 
REM  THE  ATMOSPHERE  AFTER  BLANCHARD  <1963) 

REM 

REM  FIND  THE  DROPLET  DIAMETER  AT  91.  47.  RH 
D9=P<3>*<  1+0  9/  <  1-0.  914)  )  ''<  1/3) 

W=U9(P (22) )*1.  9438 
IF  W<4.  5  THEN  20700 
C2=-13.  3+8.  966*L0G  <  W ) 

C2=C2/ 1000 

C3=14.  301-0.  386*W 

C3=-C3/10 

C4=l .  764+0.  1 91  *W 

F=C2*EXP(C3*L0G(C4/<0.  5*D9) )~2) 

F=1737.  178*F/D9 
P  ( 1 5 )  =F 
GO  TO  20705 
P< 15>=0 
GOSUB  20740 

F*= " ©AEROSOL /PROGR AM/ TASKS " 

GO  TO  30060 

REM  ROUTINE  TO  FIND  THE  MAX  TIME  STEP  ALLOWED 
P < 19 ) =60  MIN  1/(1.  0E-5+K(l>) 

FOR  1=1  TO  23 
Z1=Z<I+1)-Z(I) 

Z2=Z 1^2 

K1  =  (K( I >+K( 1  +  1 ) )/2+l.  0E-3 
T9=Z2/(2*K1) 

P<19)=P<19)  MIN  T9 
NEXT  I 
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20820  P< 19)=P< 19) -0  1 
20830  RETURN 
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20000  REn 

20010  REM  J - ! 

20020  REM  !  ©AEROSOL /PROGRAM /TASK 5  l 

20030  REM  {  S  GATHMAN  NRL  CODE  4327  J 

20040  REM  i - ! 

20050  REM 
20060  SET  KEY 
20070  REM 

20080  REM  GO  TO  PARTIAL  DIFFERENTIAL  ANALYSER 
20090  GOSUB  20660 
20100  REM 
20110  REM 

20120  REM  GO  TO  CLOCK  UPDATE  SUBROUTINE 
20130  GOSUB  20940 
20140  REM 

20150  REM  GO  TO  STORAGE  ROUTINE 
20160  GOSUB  21170 
20170  REM 

20180  REM  PARAMETER  UP  DATE  ROUTINE 
20190  IF  P  < 1 7 ) =0  THEN  20240 
20200  P  < 17 ) =0 

20210  F**M®AERaS0L/PR0GRAM/TASK2” 

20220  SET  NOKEY 
20230  GO  TO  30060 

20240  GO  TO  P(20)  OF  20090.20310.20360,20490,20580 
20250  END 
20260  REM 

20270  REM************************************************************* 
20280  REM  ******************KE YBOARD  SERVICE  ROUTINES***************** 
20290  REM  *********************************************************** 
20300  REM  PRINT  MODEL  TIME  (SEC) 

20310  PRINT  P ( 24 ) 

20320  P ( 20 ) =1 
20330  GO  TO  20090 
20340  REM 
20350  REM  PLOT  N1 
20360  PAGE 
20370  N9*-l.  0E+24 
20380  FOR  1*1  TO  24 
20390  N9*N9  MAX  N1 < I) 

20400  NEXT  I 

20410  WINDOW  O#  N9#  0>  Z  <  20 ) 

20420  AXIS 

20430  MOVE  N1<1>, Z(l) 

20440  DRAW  Nl,  Z 
20450  P ( 20 ) = 1 
20460  GO  TO  20090 
20470  REM 
20480  REM  PRINT  Nl 
20490  PAGE 

20500  PRINT  "ALTIM)",  "#/M~3",  MK<Z>” 

20510  FOR  1*24  TO  1  STEP  -1 
20520  PRINT  Z < I ) , Nl ( I ) , K ( I ) 

20530  NEXT  I 
20540  P ( 20 ) * 1 
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20550  GO  TO  20090 
20560  REM 

20570  REM  SAVE  THE  WORLD  ROUTINE 
20580  CALL  "REWIND",  2 
20590  P ( 20 ) = 1 

20600  WRITE  #2  P ,  N 1 ,  K,  H9,  U9,  2  ,  A i  ,  C 1 ,  D 1 ,  SO,  R8,  V 
20610  INIT 
20620  END 
20660  REM 

20670  REM========PARTIAL  DIFFERENTIAL  EQUATION  ANALYSER 

20680  REM 
20690  T1=P< 19) 

20700  H=P ( 22 ) 

20710  S6=  Z ( 2 ) — Z ( 1 ) 

20720  S7=(Z(3)-2U  )  )/2 
20730  Sl=-(01 (2)-0l < 1 ) )/S7 
20740  N 1 ( 1 ) =01 ( 1 ) +T1* < P ( 1 5 ) — K ( 2 ) *S1 ) /S6 
20750  N1 < 1 >=N1 < 1 )+V( 1 )*T1*(Q1 <2)-01 ( 1 ) ) /S6 
20760  FOR  1=2  TO  22 
20770  S7=Z  < I +1 ) -Z  < I ) 

20780  S8= (Z( 1+2 )-Z  < I ) )/2 
20790  S9=<Z< 1+1 >-Z< 1-1 ) >/2 
20800  Sl=-<01 < 1+1 )-01 < I ) ) 

20810  Rl=-(01 < I )-01 < 1-1 ) ) 

20820  N1 ( I ) =01 < I ) +T 1  *  <  K ( I )*R1/S9-K< 1  +  1 )*S1/S8)/S7 
20830  N1 ( I )-Nl  < I )-V< 1 )*T1*S1/S7 
20840  NEXT  I 
20850  FOR  1=1  TO  24 
20860  IF  N 1 < I ) =>0  THEN  20880 
20870  N 1 < I ) =0 
20880  NEXT  I 
20890  01 =N1 
20900  RETURN 
20910  REM 

20920  REM= ========= =====CLOCK  UPDATE  SUBROUT  I NE  =  ===  =  =  =  === 

20930  REM 

20940  P  <  24 ) =P ( 24 )  +P  ( 19) 

20950  P<23)=P(23)+P< 19) 

20960  IF  P  <  23 ) <600  THEN  21010 
20970  P  <  23 ) =0 
20980  P ( 21 >  =  1 
20990  P(18)=P<18)+1 
21010  P< 16)=P( 16)+P< 19) 

21020  IF  P ( 1 6 ) <3600  THEN  21110 
21025  P ( 22 ) =P ( 22 ) + 1 
21030  P ( 1 6 ) =0 
21040  P< 1 7 ) = 1 

21110  IF  P  < 1 8 ) < 1 44  THEN  21130 
21 120  P  <  20 ) =5 
21130  RETURN 
21140  REM 

21150  REM=====  =  =  =  =  ===  =  =  5T0R AGE  SUBROUT  I NE  =  ===  =  =  =  =  =  =  =====  = 

21160  REM 

21170  IF  P ( 2 1 ) =0  THEN  21200 
21 180  P  <  21 ) =0 
21190  WRITE  # 1 #  P ( 1 8 )  N1 
21200  RETURN 
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